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The extragalactic background light at far-infrared wavelengthi^^ originates from optically- 
faint, dusty, star-forming galaxies in the universe with star-formation rates at the level of 
a few hundred solar masses per year"*. Due to the relatively poor spatial resolution of far- 
infrared telescopes^ ^5 the faint sub-millimetre galaxies are challenging to study individually. 
Instead, their average properties can be studied using statistics such as the angular power 
spectrum of the background intensity variationiP'^. A previous attempt^^ at measuring this 
power spectrum resulted in the suggestion that the clustering amplitude is below the level 
computed with a simple ansatz based on a halo modeP^. Here we report a clear detection of 
the excess clustering over the linear prediction at arcminute angular scales in the power spec- 
trum of brightness fluctuations at 250, 350, and 500 /xm. From this excess, we find that sub- 
millimetre galaxies are located in dark matter halos with a minimum mass of log[Mniin/M0] 
= ll.Slg^ at 350 /xm. This minimum dark matter halo mass corresponds to the most ef- 
ficient mass scale for star formation in the universeP!, and is lower than that predicted by 
semi-analytical models for galaxy formatioii^. 

Despite recent successes in attributing most of the extragalactic background light at sub- 
millimetre wavelengths to known galaxy populations through stacking analyse^^^H^, we have not 
individually detected the faint galaxies that are responsible for more than 85% of the total extra- 
galactic intensity at these wavelengthd^. The faint star-forming galaxies are expected to trace the 
large-scale structure of the Universe, especially in models where galaxy formation and evolution 
is closely connected to dark matter halos. While not individually detected in low resolution ob- 
servations, the clustering of galaxies is expected to leave a distinct signature in the total intensity 
variations at sub-millimetre wavelengths. The amplitude of the power spectrum of intensity vari- 
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ations as a function of the angular scale provides details on the redshift distribution and the dark 
matter halo mass scale of dusty, star-forming galaxies in the universeP. 

For this analysis, we used data from the Herschel Multi-tiered Extra-galactic survey (HerMEI^^, 
taken with the Spectral and Photometric Imaging Receiver (SPIRE^^) onboard the Herschel Space 
Observator)^, during the Science Demonstration Phase (SDP) of Herschel. The data are com- 
posed of a wide 218' by 218' area in the Lockman Hole complemented by a narrow, but very deep 
(30 repeated scans), map of the Great Observatories Origins Deep Survey (GOODS) North field 
covering 30' by 30'. These fields have been very well studied at other wavelengths and they are 
known to have a low Galactic dust density, making it easier to distinguish the extragalactic com- 
ponent we wish to study. The observing time to complete each of the two fields was about 13.5 
hours, observing simultaneously at 250, 350, and 500 /xm. 

To limit the influence of a few bright galaxies on the measurement of the power spectrum, 
we removed galaxies brighter than 50 mJy in all three passbands by masking pixels in our maps 
with values larger than 50 mJy/beam as well as the neighboring pixels. We use the cross-power 
spectrum of two sub-maps as our estimate of the sky power spectrum to remove the contribution 
from the instrumental noise and alleviate potential systematic effects. We correct the raw cross- 
power spectra for the effects of the angular response function of the instrument and the transfer 
function of the map-making process. The angular response is established from a different set of 
SPIRE observations targetting Neptune, a strong point-like source for SPIRE, and involving a fine 
sampling of the beam with a total of 700 scan^. The effects of filtering of the time-ordered 
data and of the map pixelization are captured with a large set of sky simulations. To estimate our 
uncertainties, we propagate the errors from the beam measurement, while the simulations provide 
us with the instrumental and sky variance. The quadratic sum of these errors constitutes our error 
estimate. 

The measured angular power spectrum (c.f. Figure la) contains contributions from spatial 
variations in Galactic dust clouds, or cirrus, brightness at large angular scales, the clustering of 
galaxies at intermediate angular scales, and a white-noise component at small angular scales arising 
from the Poisson behavior of the faint galaxies^ ^. The cirrus signal in our Lockman-SWIRE field 
is taken from existing measurements in the same field with IRAS 100 /xm and MIPS^^ and extend 
this spectrum from 100 /xm to SPIRE wavelengths using the spectral dependence of a Galactic 
dust modeP^. We remove this cirrus power spectrum from our measurements and account for the 
uncertainty of the cirrus power spectrum by adding its error in quadrature to errors of our power 
spectrum points. 

The Poisson behaviour of sources lead to an additional term to the angular power spectrum 
that is scale independent. The clustering component we measure is the excess of clustered back- 
ground fluctuations above this shot-noise level. As the confusion noise is at the level of 6 mJy 
at SPIRE wavelenghtP, with fluctuations in the brightness of the background we are probing the 
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clustering of faint galaxies with fluxes at the level of a few mJy at 350 /xm. To extract astrophys- 
ical information on faint galaxies from the clustering power spectrum we make use of the halo 
modep21. This phenomenological approach connects the spatial distribution of galaxies in the uni- 
verse to that of the dark matter. To model sub-millimetre galaxies in dark matter halos, we take a 
halo occupation distribution describing the number of galaxies as a function of the halo mass M of 
A^gai(M) = 1 + {M/MiY when M > M^in above the minimum mass scale Mmin, Mi is the mass 
scale at which more than one galaxy is present in a single dark matter halo, taken to be between 
10 to 25 times M^in, and a is the power-law scaling of satellite galaxies with halo mass. The halo 
model involves two parts with clustering of galaxies within halos, the one-halo term, and cluster- 
ing of galaxies between halos, the two-halo term. While with the two-halo term alone parameters 
related to the occupation number are degenrate with each other and the bias factor or the number 
density of galaxies, with clustering in the one-halo part of the power spectrum aslo included the 
parameter degeneracies are broken and M^in can be determined more accurately^. 

At scales of a few arcminutes and above we measure a clustering excess arising from the 
one-halo term and above the two-halo term tracing the linear density field power spectrum scaled 
by galaxy bias (c.f. Figure lb). The one-halo term arises when more than one far-infrared galaxy 
occupies the same halo. While a hint of this one-halo term was previously seen in the clustering of 
the bright (> 30 mJy) sub-millimetre galaxie^, evidence for such clustering was not present for 
bright galaxies in a different Herschel dataset^"^. To describe the power spectrum of the intensity 
fluctuations we also need a prescription for the redshift evolution of the source intensity. While 
models exist in the literaturePSI^^, our data are of sufficient quality that we can directly constrain the 
redshift evolution of the source intensity from our measurements and we constrain its value in four 
redshift bins in the range to 4.0. The halo model parameters, the source intensity parameters, and 
the shot-noise amplitude are jointly estimated with Markov-Chain Monte-Carlo fits to the power 
spectrum measurements. We impose an additional prior on our parameter estimates that the redshift 
integrated source intensity from our model fits, including the fractional contribution from bright 
sources that we have masked, be within the 68% range of the known background light intensity at 
each of the three wavebands^. We combine the estimates of the source intensity evolution at three 
wavebands of 250, 350 and 500 /xm and in four redshift bins to measure the bolometric luminosity 
density between 8 and 1100 /xm as a function of the redshift (c.f. Figure 2). We find that the 
far-infrared luminosity density continues to be significant out to a redshift of 4 and is at least a 
factor of 10 larger than the luminosity density of individually detected sub-mm sources with flux 
densities above 30 mJy alonePl. 

Using the halo model fits, we estimate the mininum dark matter mass scale for dusty star- 
forming galaxies at the peak of the star formation history of the universe to be log^Q Mmin/M© = 
ll-5ta2 350 jim with a bias factor for the galaxies of 2.4to:2- The minimum halo masses 
log^Q Mmin/M© at 250 and 500 /xm are ILltg^e H-Sioi' respectively. The corresponding bias 
factors for the galaxies are 2.0to:i and 2.8to:5 at 250 and 500 /xm, respectively. The differences 
in the minimum halo masses and the bias factors between the three wavelengths are likely due 
a combination of effects including overall calibration uncertainties, the fact that at longer wave- 
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lengths we may be probing colder dust than at shorter wavelengths, and due to differences in the 
prior assumption on the total background intensity. In future, numerical models on the distribution 
of sub-millimetre galaxies will become useful to properly understand some of these subtle differ- 
ences. Averaging over the three wavelengths the minimum halo mass for sub-millimetre galaxies is 
at the level of 3 x 10^^ M©, with an overall statistical uncertainty of roughly ±0.4 in the logarithm 
of the minimum halo mass. 

Based on a variety of observed scaling relations such as the one between stellar mass and 
circular velocity, the dark matter halo mass scale for efficient star-formation has been indirectly 
inferred to be about 10^^ M(^. As the sub-millimetre galaxies are the most active star-forming 
galaxies in the universe, it is likely that the minimum halo mass scale for such galaxies that we 
determine from brightness fluctuations corresponds to the preferred mass scale of active star- 
formation in the universe. In dark matter halos below this mass, star-formation is expected to 
be inefficient due to photo-ionization feedback^^. The underlying astrophysics needed to explain 
the numerical value we find are still missing in galaxy formation theories, since existing semi- 
analytical models predict a mass scale for faint sub-millimetre galaxies that are roughly ten times 
largei^. We provide the strongest evidence for a minimum mass scale for active star-formaing 
galaxies in the universe by stuyding the background intensity variations generated by those galax- 
ies in our sky maps. Our direct estimate of the minium dark matter halo mass provides a critical 
value needed to improve theoretical models of sub-millimeter galaxies and the overall galaxy for- 
mation and evolution. 
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Figure 1: The two-dimensional power spectrum of the Herschel map. (a): The total power 
spectrum P{k) at 350 /xm as a function of the wavenumber k in arcmin"^. The error bars are 
68% confidence level uncertainties. The shaded region is the cirrus signal in our Lockman-SWIRE 
field. To describe the total power spectrum, we take a power-law with P{k) = A{k/ki)^ + Psn 
where ki is fixed at 0.1 arcmin"^ and Psn is the shot-noise amplitude. At 350 /xm, we find A = 
(5.79 ± 0.26) X 10^ Jy2/sr and n = -1.28 ± 0.07, shown by the red-line. The light blue line 
is the best-fit shot-noise amplitude with a value of 4600 ± 70 Jy^/sr, in agreement with the value 
of 4500 ± 220 Jy^/sr predicted by best determined source countP^. The shot-noise errors include 
the 15% uncertainty in the SPIRE absolute flux calibratiorPH. in green we show the total power 
spectrum combined with the mean estimate of the cirrus signal, (b): The angular clustering power 
spectrum at 350 /xm as a function of the wavenumber with errors showing the 68% confidence 
level uncertainties. The best-fit shot-noise value (dashed blue line) has been removed from the 
data and its uncertainty added to the overall error in quadrature. The green lines show the best- 
fit halo model with a reduced chi-square value of 1.02. The dot-dashed line show the two-halo 
term and the dashed line show the one-halo term that is responsible for the clustering at small 
angular scales. Data in the lowest wavenumber bins contaminated by the Galactic cirrus have 
been omitted. For comparison we also show a previous measurement of the power spectrum of 
brightness fluctuations at 350 /xm with BLAST^. The results related to 250 and 500 /xm are 
summarized in the Supplementary Information. 



Figure 2: Far-infrared bolometric luminosity density (between 8 and 1110 /xm) and star for- 
mation rate as a function of redshift. The far-infrared bolometric luminosity density was esti- 
mated using the values of the source intensity in four redshift bins between and 4 derived from 
best-fit model fits to each of the power spectra from the three wavebands and with a prior selection 
of models with a > 1 for the power-law slope of the occupation number. We have assumed a mod- 
ified black-body spectrum for the spectral energy distribution of the spatially unresolved sources, 
with an emissivity index of 1.5 and a dust temperature of 28±8 K^. We propagate the uncertainty 
on the dust temperature as an additional error when computing errors on the luminosity density 
in each of the redshift bins. The vertical bars show the 68% confidence level errors propagated 
from errors in the source intensity evolution estimates and from the temperature prior, while the 
horizontal bars indicated the redshift range of each of the bins we use for model fits. The orange 
and blue shaded areas represent the model from Lagache et al.^^ and Valiante et al.^^ with the same 
flux cut (S> 50 mJy) and the same temperature prior as our data. Our measurements are consistent 
with both these models. The redshift evolution of the far-infrared luminosity is a measure of the 
star formation history of the universeP^l and converting our estimate, we find a star formation rate 
of < 0.01, 0.05 ± 0.03, 0.04 ± 0.03, and 0.06 ± 0.04 M^/yr/Mpc^ in each of the redshift bins of 
z<l, l<z<2, 2<2:<3, and 3 < z < 4, respectively. 
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In these Supplemental Notes to our main paper, we outline key details related to how we 
estimated the angular power spectrum of Herschel-SPIRE data and its interpretation. 



The data presented in the main paper are publicly available from the ESA/Herschel Science 
Archive (http://herschel.esac.esa.int) under the observational identifications 1342186108, 1342186109, 
and 1342185536. Derived products by the HerMES collaboration, such as source catalogs, will be 
released through the HeDaM Database (|http : //hedam . oamp . f r/HerMES]) . 



Data Analysis and Map making: The data were taken in the standard Scan-Map AOT for which 
the scanning speed is 30"/s. Calibrated time-ordered data were created using HlPriSl develop- 
ment version 2.0.905, with a fix applied to the astrometry (included in more recent versions of 
the pipeline), with newer calibration files (SPIRE Beam Steering Mirror calibration version 2, flux 
conversion version 2.3 and temperature drift correction version 2.3.2) and with a median slope 
subtracted from each timeline. We removed a few percent of the data samples which were contam- 
inated by cosmic-rays or instrumental effects and flagged by the initial pipeline. 

We convert time-ordered data to a map on the sky through an iterative baseline removal and 
an iterative calculation of detector weights. We give a summary of the map-making method here. 
Full details of this approach is available elsewhereP"^. 

Given sky brightness I {9), the signal for a detector in a scan s with a time sample j can be 
written as 

Sdsj = H^dsj) + Pds + ^dsj -> (1) 

where P^^ is an n-th order polynomial baseline offset for detector d and scan s, and N^sj is the 
instrumental noise. The parameters of P^^ are solved with an iterative solution to the best-intensity 
of the sky at each step. 



At each iteration z, we minimize the variance of the residual V^^- based on the previous map 
I'~^{Odsj) such that 

+ • (2) 



The i-th estimate of the sky is performed via the weighted mean of all the samples that fall into a 
given pixel: 



Yjdsj ^ds i^^dsj Pds 



^dsj ^ds 

The weight associated with each iterative estimate is simply the inverse variance of the residual 

A^tot 



(3) 



w 



ds 



(4) 
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The algorithm also allows us to create a noise map by propagating detector noise as estimated 
by the variance of the residuals. The maps used here make use of a first-order polynomial with 
n = 1 and 20 iterations. The gain offsets were computed from the 1st iteration while the weights 
are fixed to 1.0 for the 10 iterations and are calculated from the data starting at iteration 11, to 
improve stability of the algorithm. The same procedure is repeated for both Lockman-SWIRE and 
GOODS-N maps and also iterative maps of Neptune that we use for beam measurements. The 
maps at 250 /xm are shown in Figure SI. 

The absolute astrometry of the maps was corrected by stacking at the positions of Spitzer 
Multi-Band Imaging Photometer (MIPS) 24 /xm and radio sources, finding reasonably consistent 
results between the two within 0.5" . We have made an overall correction to the absolute astrometry 
of the order of a few arcseconds, though such small angular scale corrections do not impact results 
we present here focusing between 0.5 arcminutes to 100 arcminutes. 

The maps were made with pixels of size 6, 8.3, and 12 arcseconds at 250, 350, 500 /xm, 
respectively, corresponding to one third of the full- width at half-maximum (FWHM) of the SPIRE 
beam profiled in the three passbands. 

Raw Power Spectra To compute the power spectrum in our map, we use fast Fourier transforms; 
however, we need to take into account the missing and unwanted pixels. Due to our scanning strat- 
egy and some corrupted data, a small fraction of pixels are not defined on the map that we use to 
define our Fourier transform basis. Furthermore, we wish to remove the brightest galaxies in order 
to reduce the shot-noise term in the power spectrum, since the shot-noise term is weighted towards 
bright galaxies and larger shot-noise degrades the ability to extract the clustering component of 
the power spectrum. In our analysis we applied a flux cut of 50 mJy/beam, removing 0.9, 0.7, 
1.2% of the pixels at 250, 350, 500 /xm, respectively, in the Lockman-SWIRE field and 0.5, 0.2, 
0.2% in the GOODS-N field. We used the same flux cut at all frequencies for simplicity, this 50 
mJy/beam allows to remove all the bright sources while retaining most of the pixels in each map. 
The remaining number of pixels used for the fluctuation study is 5.4 x 10^ , 2.9 x 10^, and 1.4 x 10^ 
at 250, 350, 500 /xm, respectively, in the Lockman-SWIRE field and L9 x 10^ LO x 10^ 4.7 x 10^ 
in the GOODS-N field. 

The raw power spectra are summarized in Figure S2. Here, we show the auto spectra in the 
total map as well as the cross spectrum with maps made with half of the time-ordered data in each 
map. The difference of the two provides us with an estimate of the instrumental noise. At small 
physical scales (large k values) the noise is almost white such that P{k) is a constant value. We fit 
a model of the form 



and determine the knee scale of the noise, fco, to be at about 0.15 arcmin ^ at each of 250, 350, and 
500 /xm. 




(5) 
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Figure S 2: Raw P{k) of the Lockman-SWIRE (left) and GOODS-N (right) fields. The panels 
show the three wavebands with 250 jim (panels a and d), 350 fim (b and e), and 500 /xm (c and 
f) from top to bottom on each side, respectively. The blue lines show the auto-power spectra 
computed from maps using all the data. This power spectrum is a combination of sky signal and 
instrumental noise. We estimate the sky signal through the cross-spectrum (green lines) of two 
maps after dividing the data into two halves separated in time. The difference of these two spectra 
represents an estimate of the instrumental noise power spectrum (red line). 
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Corrected Power Spectra: The final power spectrum we show in this paper is corrected for a 
combination of effects described by 

P{k') = B{k)T{k)M,.,P{k) , (6) 

where P{k') is the observed power spectrum from data in the presence of mask, B{k) is the beam 
function, and the map making transfer function is T{k). P{k) is the true sky power spectrum and 
is determined by inverting the above equation. 

In the above, M^k' is the mode couphng matrix associated with the mask. This can be 
expressed analytically in the flat sky approximation^^ as 

Mkk^ = E E - k')\^/N{9k) , (7) 

Ok Oy 

where w{k) is the Fourier transform of the mask. Figure S|3] shows the mask we used and the 
corresponding matrix M^^f for each of Lockman-SWIRE and GOODS-N fields at 250 /xm with 
sources above 50 mJy and spurious data removed. 

We measure the beam function B{k) through observations of Neptune, involving a total of 
700 scans. Figure S|4] shows the Neptune maps made at 250, 350, and 500 /xm. The Neptune 
data are analyzed in the same manner as the Lockman-SWIRE and GOODS-N data using the 
same iterative map maker. In addition we also employ a naive map maker available as part of 
HIPE^^ When making these maps, we account for the relative motion of Neptune relative to the 
background sky and make maps that correct for Neptune's varying position during the observations. 
This results in a map where extragalactic sources are smeared. Neptune, however, is several orders 
of magnitude brighter and our beam measurements primarily focus on the central region. Figure S|5] 
summarizes the results related to B{k) for each of the three SPIRE wavelengths. In the same figure 
(bottom panels), we also compare the beam measured from Neptune to the beam described by a 
Gaussian with a FWHM of 18, 25, and 36'' at 250, 350, and 500 /xm. The amplitude of B{k) is 
thereafter interpolated in the k modes at which we compute our fluctuation power spectra. 

The uncertainty in the beam function B{k) is determined by computing the standard devi- 
ation of the different B{k) estimates, using the measurements on the iterative and naive map and 
several different interpolation schemes. The beam uncertainty computed this manner is slightly 
larger than the difference in the beams between two different observations of Neptune, one involv- 
ing the fine scans we primarily use here and an older coarse set of Neptune scans, but with maps 
made using the same map-maker. Figure S|6] shows the overall uncertainties in the beam (solid 
lines) as well as the uncertainties coming from the difference between the naive and iterative map 
reconstructions (dashed lines). Figure S|7] compares the beam function and the power spectrum 
P{k) at 250 /xm from the Lockman-SWIRE field showing that features in the power spectrum are 
not related to features in the beam function B{k). 

To measure the transfer function T{k) associated with the map maker, we realize 100 simu- 
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lations of a first estimate of our beam-convolved power spectrum and pass it through the iterative 
map-making pipehne used to reduce our real data. We then compute the average of the ratio be- 
tween the estimated spectrum and the input spectrum with simulated maps masked exactly as in the 
real sky maps. This function is the transfer function T{k) of the map-making pipeline associated 
with median filter and other filtering (Figure S©. We divide the estimated power spectrum by this 
transfer function to remove the map making pipeline processing effects. 

Total error budget: The total error budget in our clustering plots is composed of three contribu- 
tions involving the uncertainty of the beam, uncertainty in the shot-noise determination, and the 
instrumental noise error. The latter is computed from simulations while the beam error comes 
from the differences in our estimates of the beam. Shot-noise error results from direct fits to the 
measured data points. Figure S|9] summarizes the error budget as a function of wave number and 
also compares the instrumental noise error to an analytical formula for its expectatioiP^. 

While in Figure 1 of the main paper we showed P{k) at 350 /xm, in Figure SfTOlwe show the 
same at 250 and 500 /xm. 

Jack-knife tests: The results shown in this study were obtained by dividing the data into two 
equal and consecutive halves and by taking the cross-power spectrum of the resultant maps. We 
can make the same cross-correlation power spectrum measurements and repeat the whole process 
with maps made by dividing data into several other combinations. To be specific, we divide the 
data into four pieces, each filling approximately our total field, and use these to measure the cross- 
power spectrum for two other combinations namely [(1 + 3) x (2 + 4)] and [(1 + 4) x (2 + 3)], 
where 1 to 4 are four equal subdivisions of data in time. 

Figure SfTT] summarizes our results, showing that within the uncertainties we recover similar 
power spectra. Given the observing strategy, the (1 + 3) and (2 + 4) maps are each made with 
parallel scans, but roughly perpendicular to each other. The fact that we do not see a statistically 
significant difference shows that the beam ellipticity is not an important systematic concern in this 
study. 

Null tests: In addition to the jack-knife tests with a variety of sub-maps with data divided to four 
intervals and all leading to a measurement of the sky signal, we also perform several null tests 
using data combinations that remove the sky signal. In this case, instead of computing the cross- 
power spectra of the sum maps of data combinations, we make use of the sub-maps made by taking 
the differences of data combinations, again data divided to four sub-intervals as the case of signal 
measurement. As an example, in Figure SfT2lwe show the cross-power spectrum computed at 250, 
350, and 500 fim with the (1 — 2) map cross-correlated against the (3 — 4) map. For reference, we 
also show the default power spectrum computed with [(1 + 2) x (3 + 4)]. 
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Figure S 3: Mask and the model coupling matrix. Top-left : Mask used for the Lockman- 
SWIRE field. The galaxies brighter than 50 mJy masked as well as some corrupted scans. Top- 
right: Coupling matrix M^k' (log scale) computed for this Lockman-SWIRE mask. Bottom-left 
and bottom-right figures are the same things for the GOODS-N field. 



18 





0123012301234 

Figure S 4: Maps used for the beam function measurements. From left to right, Neptune beam 
maps (log scale normalized to the maximum) at 250, 350, and 500 /xm, respectively. The x- and y- 
axes are labeled in arcminutes. The maps were made with the iterative map-maker and observations 
involve a total of 700 scans of Neptune. 
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Figure S 5: Point spread function of SPIRE Instrument. Top : Point Spread Function Fourier 
space kernels for 250, 350 and 500 /xm (from left to right). The black diamonds were estimated 
using the Neptune "naive" map, the blue triangles with the Neptune "iterative" map and the red 
dashed line is the Gaussian (FWHM of 18, 25 and 36 arcseconds). The vertical green dashed dotted 
lines represent the maximum k out to which the data are used in this analysis. Bottom: Ratio of 
the beam kernel measured on Neptune to the Gaussian beam approximation. 
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Figure S 6: Accuracy of the beam measurement. The beam uncertainty relative to the mean 
beam function used for this analysis. The total uncertainty (solid lines) is the standard deviation 
of the different B{k) estimates using the measurements on the iterative and naive map and several 
interpolation methods. The dashed lines are half of the difference between one of the naive map 
B{k) estimates and one of the iterative map B{k) estimates. The vertical lines mark the maximum 
k value out to which we make use of the power spectrum measurements at each of 250, 350, and 
500 fim. 



Band 


A (Jy2/sr) 


n 


PsN (Jy'/sr) 


xVd.o.l 


250 jum 
350 fim 
500 


(7.64 ± 0.55) X 10^ 
(5.79 ± 0.26) xlO^ 
(2.67±0.13)xl03 


-1.20 ±0.09 
-1.28 ±0.07 
-1.16 ±0.09 


57981^^2 
43731?^ 
1700 ± 80 


0.93 
1.03 
1.2 



Table S 1 : Power-law best fit values at fci = 0.1 arcmin 



Notes: To describe the power spectrum, we take a power-law with P{k) = A{k/kiY + ^sn 
where ki is fixed at 0.1 arcmin"^ and Psn is the shot-noise amplitude, assuming a power- 
law fit to the data. The errors are 68% confidence level uncertainties determined from the 
MCMC fits. 
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Figure S 7: Power spectrum relative to the beam function. Comparison of the P{k) estimate 
and the beam transfer function shape. The black diamonds represent our 250 jim P{k) estimate 
divided by the best fit power-law. The red line represents our beam transfer function divided by the 
approximate Gaussian beam. The two curves have different shapes and this difference indicates 
that the shape does not come from the beam transfer function. 
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Figure S 8: Map-making transfer function. Transfer function T{k) due to the iterative map- 
maker and filtering on the cross-spectra of our two fields. The blue, green and red triangles rep- 
resent, respectively, the transfer function at 250, 350 and 500 /xm. T{k) is essentially equal to 1 
between 0.02 and 0.4 arcmin"^. The map-maker adds about 10% power around 0.01 arcmin"^ 
in the case of Lockman-SWIRE (top panel) and around 0.05 arcmin"^ in the case of GOODS-N 
(bottom panel) and reduces the power on small scales mostly by averaging the data into pixels 
(light blue dotted dashed lines). The vertical lines mark the maximum k out to which we make use 
of the power spectrum estimates for shot-noise and clustering measurements. 
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Figure S 9: The power spectrum uncertainties. Error budget at 250, 350, and 500 /xm. We show 
the error separated into the beam uncertainty (blue hues), the shot-noise determination (green 
hues), and the simulations (sky and instrumental variance, red lines). The simulation uncertainty 
is compared to an analytical noise estimate (black lines). 
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Cirrus power spectrum The cirrus signal in our Lockman-SWIRE field is taken from existing 
measurements in the same field with IRAS 100 /xm and MIPI^ with a power spectrum, P{k), of the 
formPcirrus(fc) = Pq (fc/^o)^ at 160 /xm with Pq = (2.98 ± 0.66) x lOMy^/sr and /3 = -2.89±0.22 
when ko = 0.01 arcmin"^. We and extend this spectrum from 100 /xm to SPIRE wavelengths using 
the spectral dependence of a Galactic dust modeP^. 
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log[M^jMQ] 


a 




PsN (Jy'/sr) 


x'/d.o.f. 


250 fim 




-■-•0-0.2 


2.01°:? 


6100 ± 120 


0.76 


350 




-L-^-O.T 


2^41S:^ 


4600 ± 70 


1.02 


500 fxm 


11.81°:! 


1 8+0-1 
-•-•0-0.7 


2.8l°:t 


1800 ± 80 


1.44 



Table S 2: Halo model best fit values from the measured power spectra at the three wave- 
bands. 

We tabulate the best-fit values with 68% confidence level errors for halo occupation num- 
ber used to interpret the power spectrum measurements. The average galaxy bias factor 
is (6)^. PsN is the amplitude of shot-noise fluctuations, also jointly determined from the 
power spectra as part of our model fitting process. The errors of the shot-noise amplitudes 
PsN include an extra error corresponding to the uncertainty of the absolute flux calibra- 
tion scale at the three SPIRE wavebands of 15°/cP. The chi-square values of the best-fit 
model, per degree-of-freedom, are also tabulated. We do not tabulate the values of Mi 
as it remains unconstrained within the prior of Mi/M^m taken to be between 10 and 25. 
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Interpretation Model In Table SI we summarize results related to power-law model fits with 
P{k) = A{k/ki)^ + PsN. where /ci = 0.1 arcmin"^. We also make use of the halo model approach 
to model-fit our clustering measurements, making use of the halo occupation distribution (HODj^. 

The number of pairs of galaxies inside a given halo depends on the variance of the HOD, 
cr^(M, z) = ( A^gai(A^gai — 1)) whilc the number of pairs of galaxies in different halos is simply 
given by the square of the mean halo occupation with N{M, z) = ( A^gai), where A^gai is the total 
number of galaxies in a halo and we further assume that one galaxy occupies the center of the halo, 
the others being considered as satellite galaxies, so that A/gai = A^cen + ^sat- Central and satellite 
galaxies are assumed to have different HODs; in fact the mean number of central galaxies in a 
given halo is a simple step function, so that A^cen = 1 above a given mass M^in, and A^cen = 
otherwise. The HOD of satellite galaxies is taken to be a power law of the halo masd^ 

^sat = hr:r ; (8) 



Ml 

here Mi is a normalization factor that represents the mass scale at which a single halo hosts on 
average one satellite galaxy in addition to the central galaxy. 

The power spectrum of galaxies is then parameterized as the sum of two different contri- 
butions: the 1-halo term, which describes the clustering on small scales and is related to pairs of 
galaxies within the same halo and the 2-halo term, responsible for the large scale clustering and 
related to pairs of galaxies in different halos: 

P(fc,z) = Pih(fc,z) + P2h(fc,z). (9) 

The two terms are then written as: 

Nl,{M)uly,{k,z\M)]dM/nl,,{z), 
P2h(k,z) = PBM(k,z)x 

rf^halo 



JdM^{z)N,^,{M,z) 



2 



b{M,z)uBM{k,z\M)dM\ 

(10) 

Here Prtuik, z) is the linear dark matter power spectrum; nhaio is the halo-mass functioiP, 6(M, z) 
is the linear bias which connects the large scale clustering of dark matter to the galaxy clustering; 
«dm(^, z\M) is the normalized dark matter halo density profile in Fourier space (as a function of 
wavenumber k and redshift z for a given value of mass M) and ngai is the mean number of galaxies 
per unit comoving volume at redshift z: 

M 



(11) 
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For the dark matter halo density function we adopt the Navarro-Frenk-WhiteP' profile truncated at 
the virial radius rvir and with a concentration parameter given by: 



^^^'^^-TnUj ^ '''' 

here is the characteristic mass scale at which the critical density required for spherical collapse 
is equal to the square root of the variance in the initial density field a{M) when extrapolated at the 
present time using linear theory such that 5sc(^) = cr(M*), where S^dz) = 1.68/ g{z), where g{z) 
is the linear theory growth function for density perturbations. 



As outlined above M^in determines both the one-halo and two-halo amplitudes, while a 
determines primarily the amplitude of the one-halo term and the overall number density of galaxies, 
which in return is connected to the amplitude of the two halo term via the halo bias factor. While 
with the two-halo term alone all these parameters are degenrate with each other and the bias factor, 
allowing only an average mass scale to be determined based on the bias factor, with one-halo term 
also included some of the degeneracies are broken and M^in and a can be determined independent 
ofbiaP. 



In the Limber approximation, the measured power spectrum of fluctuations can be expressed 
as the 2 dimensional, flux averaged projection of the three-dimensional galaxy power spectrum 

P{k, z) as: 




here dS/dz is the redshift distribution of the cumulative flux contributed by the background faint 
galaxies, dV^ is the comoving volume element, defined as dV^ = and x{z) is the comoving 

radial distance. 



In this paper we determine the redshift distribution of the intensity by binning the redshift 
range in four redshift bins between z = and z = 4 and putting constraints on dS/dz in each bin; 
the advantage of this approach is that we don't assume a particular model for dS/dz{z)\ instead, 
we let the data decide which model is more adequate. 



The method we use to constrain our parameters is based on a modified version of the pub- 
licly available Markov-Chain Monte-Carlo (MCMC) package CosmoMCpl, with a convergence 
diagnostics based on the Gelman-Rubin criterion"^^. We consider a halo model described by the 
following set of parameters: 

{dS/dz,, M^in, a, Ml, Psn} , (14) 

where, as discussed before, we bin the cumulative flux dS/ dz{z) in four redshift bins, dS/ dzi {z) {i = 
1, 2, ..4), representing the value at four redshift intervals, bin^ G [0 — 1, 1 — 2, 2 — 3, 3 — 4]. In the 
above Psn is the shot-noise amplitude which we remeasure again for the halo model fits. To obtain 
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Figure S 10: The fluctuation power spectrum and the clustering component. The total power 
spectrum P{k) (left panels) and clustering P{k) with shot-noise removed (right panels) at 250 fim 
(top) and 500 /xm (bottom), respectively. The power spectrum measurements shown are binned 
logarithmically for k > 0.03 arcmin"^, with a bin width equal to Ak/k = 0.15, and linearly for 
smaller fc, with a bin width of Ak = 4.6 x 10~^ arcmin"^. This figure is similar to Fig 1. 
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Figure S 1 1 : Accuracy of the power spectrum measurement. Ratios of the total power spectra 
P{k) estimated with different sub-maps of Lockman-SWIRE data normahzed to the default power 
spectrum shown in the main paper estimated with the (1 + 2) map cross-correlated against the 
(3 + 4) map, after the data are divided into four sequential intervals in time, labeled 1 to 4, of equal 
duration. The 1 and 3 subsets have the same scan direction and the 2 and 4 subsets have the same 
scan directions, but the two subsets are almost orthogonal. 
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Figure S 12: The null-test of the power spectrum measurement. P{k) measured (black triangle) 
on Lockman-SWIRE with the cross spectrum [(1 + 2) x (3 + 4)] at 250, 350, 500 /xm (left to right, 
top to bottom). Cross power spectrum (green diamond) of the difference [(1 — 2) x (3 — 4)]. 
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Figure S 13: The halo model parameter estimates. Bi-dimensional probability distribu- 
tion function for all the pairs associated with our halo model fits with eight parameters 
(Ml, Mjnin, a, PsN, Si, S2, S3 and 5^4) showing our constraints and the degeneracies between the 
parameters. Here we show results at 350 //m, but degeneracies of parameters related to 250 and 
500 fim model fits are similar. 
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Figure S 14: Redshift evolution of the galaxy intensities at the sub-millimetre wavelengths. 

dS/dz as a function of redshift for 4 bins in redshift and for the three wavebands of SPIRE with 
< 50 mJy. For reference, we show 2 model predictions from Lagache et alP^ and Valiante et 
al.f^for the same flux cut. We have used a prior on the occupation number slope a> 1. 

reliable model-fits to data we set a broad uniform prior on the ratio Mi/Mmin to be between 10 to 
25, consistent with numerical simulations of the halo occupation distribution which finds a value 
close to 15 for this ratio"^^. We also require that the redshift integrated source intensity be within 
the 68% confidence level ranges of the background light intensities as obtained by FIRAI^. The 
central values and errors we use are 0.85 ± 0.08, 0.65 ± 0.19 and 0.39 ± 0.10 MJy/sr at 250, 350, 
and 500 /xm, respectively. For background cosmology, we assume the concordance modeP^I. Our 
results related to the halo model fits are summarized in Table S2. 

In comparison to the shot-noise values from model fits to the power spectrum (Table SI for 
the power-law case and Table S2 for the halo model case), the shot-noise values from the best 
determined source countPgive 6900 ± 320, 4500 ± 220, and 1600 ± 100 Jy^/sr at 250, 350, and 
500 /xm, respectively. 

In Figure SfTSlwe show the two-dimensional constraints on pairs of parameters that highlight 
the degeneracies associated with this eight parameter model fit. The best-fit values and the errors 
at each of the three wavebands are show in Figure S [T4l 



32 



Additionally, we compute the far-infrared bolometric luminosity between 8 and 1100 /xm 
in each of the redshift bins from the dS/dz{z) values by modelling the flux received between a 
redshift zj^^^^ and zl^^^ in each j SPIRE bands, defined by the bandpass fj{y)\ 

dS^/dz^ = LpiK I dz^—-^^—^^ (15) 



J 270 GHz 



The temperature T is chosen to be 28 ± 8 K and the emissivity index (3 is fixed to 1 .5, we then fit for 
^FiR given the measured values and the predicted values of dS/dzi. The temperature uncertainty 
is incorporated into the LpiR error budget. We summarize our results related to LpiR as a function 
of redshift in Figure 2. LpiR is a measure of the star-formation rate withP^ 

SFR[M0yr-^] = 1.73 x IQ-^^ L[L0] . (16) 

We use this to also show the SFR implied by LpiR in Figure 2. Here we have subselected the 
models that lead to a > 1 to be consistent with the occupation numbers at other wavelenegthd^'^. 
LpiR as a function of redshift has been predicted in two analytical models of sub-millimetre galaxy 
populatioiP^I^ and we make a comparison in the same figure. 
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